In [1]:
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
from numpy.fft import fft
%matplotlib inline
In [2]:
data = sio.loadmat('sunspotData.mat')
years, sunspots = data['years'], data['sunspots']
start, end = min(years)[0], max(years)[0]
ave_sunspot = np.empty(end - start + 1)
for i, year in enumerate(range(start, end + 1)):
ave_sunspot[i] = np.mean(sunspots[years == year])
plt.plot(np.arange(start, end + 1), ave_sunspot)
plt.xlabel('t [year]')
plt.ylabel('sunspots')
Out[2]:
We can clearly see in that the number of sunspots is periodic with a period of more or less 10 years
In [3]:
plt.plot(ave_sunspot[:-1], ave_sunspot[1:])
plt.xlabel(r'sunspots$_{t - 1}$')
plt.ylabel(r'sunspots$_{t}$')
Out[3]:
In [18]:
omega= fft(sunspotsdef estimated_autocorrelation(x):
n = len(x)
variance = x.var()
x = x-x.mean()
r = N.correlate(x, x, mode = 'full')[-n:]
#assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
result = r/(variance*(N.arange(n, 0, -1)))
return result).real
In [19]:
plt.plot(np.abs(np.real_if_close(omega)))
Out[19]:
In [20]:
def estimated_autocorrelation(x):
n = len(x)
variance = x.var()
r = np.correlate(x, x, mode = 'full')[-n:]
result = r/(variance*(np.arange(n, 0, -1)))
return result
In [27]:
acorr = estimated_autocorrelation(ave_sunspot.flatten())
plt.plot(range(len(acorr)), acorr)
Out[27]:
In [ ]: